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Abstract 

In timing-based neural codes, neurons have to emit action potentials 
at precise moments in time. We use a supervised learning paradigm to 
derive a synaptic update rule that optimizes via gradient ascent the 
likelihood of postsynaptic firing at one or several desired firing times. 
We find that the optimal strategy of up- and downregulating synaptic 
efficacies can be described by a two-phase learning window similar to 
that of Spike-Timing Dependent Plasticity (STDP). If the presynaptic 
spike arrives before the desired postsynaptic spike timing, our optimal 
learning rule predicts that the synapse should become potentiated. The 
dependence of the potentiation on spike timing directly reflects the time 
course of an excitatory postsynaptic potential. The presence and am- 
plitude of depression of synaptic efficacies for reversed spike timing de- 
pends on how constraints are implemented in the optimization problem. 
Two different constraints, i.e., control of postsynaptic rates or control 
of temporal locality, are discussed. 



1 Introduction 

Experimental evidence suggests that precise timing of spikes is important in 
several brain systems. In the barn owl auditory system, for example, coinci- 
dence detecting neurons receive volleys of temporally precise spikes from both 
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ears [14] . In the electrosensory system of mormyrid electric fish, medium gan- 
glion cells receive input at precisely timed delays after electric pulse emission 
[Zj. Under the influence of a common oscillatory drive as present in the rat 
hippocampus or olfactory system, the strength of a constant stimulus is coded 
in the relative timing of neuronal action potentials [23 El E21 • In humans pre- 
cise timing of first spikes in tactile afferents encode touch signals at the finger 
tips [2HJ. Similar codes have also been suggested for rapid visual processing 
[44] . and for the rat's whisker response [88] . 

The precise timing of neuronal action potentials also plays an important 
role in Spike-Timing Dependent Plasticity (STDP). If a presynaptic spike ar- 
rives at the synapse before the postsynaptic action potential, the synapse is 
potentiated; if the timing is reversed the synapse is depressed [SU E3 13 U3 IE] ■ 
This biphasic STDP function is reminiscent of a temporal contrast or temporal 
derivative filter and suggests that STDP is sensitive to the temporal features 
of a neural code. Indeed, theoretical studies have shown that, given a bipha- 
sic STDP function, synaptic plasticity can lead to a stabilization of synaptic 
weight dynamics [2H1 S3 123 SHI EE] while the neuron remains sensitive to 
temporal structure in the input [13 E3 [23 E3 E3 120] • 

While the relative firing time of pre- and postsynaptic neurons, and hence 
temporal aspects of a neural code, play a role in STDP, it is, however, less 
clear whether STDP is useful to learn a temporal code. In order to further 
elucidate the computational function of STDP, we ask in this paper the follow- 
ing question: What is the ideal form of a STDP function in order to generate 
action potentials of the postsynaptic neuron with high temporal precision? 

2 Model 

2.1 Coding Paradigm 

In order to explain our computational paradigm, we focus on the example of 
temporal coding of human touch stimuli [SB] , but the same ideas would apply 
analogously to the other neuronal systems with temporal codes mentioned 
above [HI [23 112 E2 EB] • For a given touch stimulus, spikes in an ensemble 
of N tactile afferents occur in a precise temporal order. If the same touch 
stimulus with identical surface properties and force vector is repeated several 
times, the relative timing of action potentials is reliably reproduced whereas 
the spike timing in the same ensemble of afferents is different for other stimuli 
[26] . In our model, we assume that all input lines, labeled by the index j with 
1 < j < converge onto one or several postsynaptic neurons. We think of the 
postsynaptic neuron as a detector for a given spatio-temporal spike patterns 
in the input. The full spike pattern detection paradigm will be used in Section 
[01 As a preparation and first steps towards the full coding paradigm we will 
also consider the response of a postsynaptic neuron to a single presynaptic 
spike ( Section l3.1[) or to one given spatio-temporal firing pattern fSection l3.2[) . 
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2.2 Neuron Model 



Let us consider a neuron i which is receiving input from N presynaptic neurons. 
Let us denote the ensemble of all spikes of neuron j by Xj = {£;•, . . . , tj 3 } where 
tj denotes the time when neuron j fired its k th spike. The spatio-temporal spike 
pattern of all presynaptic neurons 1 < j < N will be denoted by boldface 

X = {Xi, . . . ,X N }. 

A presynaptic spike elicited at time tj evokes an excitatory postsynaptic 
potential (EPSP) of amplitude Wij and time course e(t — tj). For the sake of 
simplicity, we approximate the EPSP time course by a double exponential 



exp [ ) — exp [ — — 



0( 



with a membrane time constant of r m = 10 ms and a synaptic time constant 
of r s = 0.7 ms which yields an EPSP rise time of 2 ms. Here 0(s) denotes the 
Heaviside step function with 6(s) = 1 for s > and G(s) = else. We set 
eo = 1.3 mV such that a spike at a synapse with w^j = 1 evokes an EPSP with 
amplitude of approximately 1 mV. Since the EPSP amplitude is a measure of 
the strength of a synapse, we refer to also as the efficacy (or "weight") of 
the synapse between neuron j and i. 

Let us further suppose that the postsynaptic neuron i receives additional 
input I(t) that could either arise from a second group of neurons or from 
intracellular current injection. We think of the second input as a 'teaching' 
input that increases the probability that the neuron fires at or close to the 
desired firing time t des . For the sake of simplicity we model the teaching input 
as a square current pulse I(t) = I Q(t - t des + 0.5AT)9(t des + 0.5AT - t) 
of amplitude Jo and duration AT. The effect of the teaching current on the 
membrane potential is 

WteachO) = / k(s)I(t - S)ds (2) 

Jo 

with k(s) = k exp(—s/r m ) where k = 1/C is a constant that is inversely 
proportional to the capacitance of the neuronal membrane. 

In the context of the human touch paradigm discussed in section 12. 1| the 
teaching input could represent some preprocessed visual information ('object 
touched by fingers starts to slip now'), feedback from muscle activity ('strong 
counterforce applied now'), cross-talk from other detector neurons in the same 
population ('your colleagues are active now'), or unspecific modulatory input 
due to arousal or reward ('be aware - something interesting happening now'). 

The membrane potential of the postsynaptic neuron i (Spike Response 
Model citeGerstner02) is influenced by the EPSPs evoked by all afferent spikes 
of stimulus x, the 'teaching' signal and the refractory effects generated by 
spikes t{ of the postsynaptic neuron 

N 

Wi(t|x, y\) = u vcst + w ij e ( f ~ tf) + Vit- t{) + iWh(*) (3) 

j=1 i^i i^vi 
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where u rest = —70 mV is the resting potential, y\ — {tj, tf , . . . , tf < t} is the 
set of postsynaptic spikes that occurred before t and tf always denotes the last 
postsynaptic spike before t. On the right-hand side of Eq. i](s) denotes 
the spike-afterpotential generated by an action potential. We take 



where r/o < is a reset parameter that describes how much the voltage is reset 
after each spike; for the relation to for integrate- and-fire neurons see [21]. The 
spikes themselves are not modeled explicitly but reduced to formal firing times. 
Unless specified otherwise, we take r] = —5 mV. 

In a deterministic version of the model, output spikes would be generated 
whenever the membrane potential Ui reaches a threshold d. In order to account 
for intrinsic noise, and also for a small amount of 'synaptic noise' generated by 
stochastic spike arrival from additional excitatory and inhibitory presynaptic 
neurons which are not modeled explicitly we replace the strict threshold by 
a stochastic one. More precisely we we adopt the following procedure pT] . 
Action potentials of the postsynaptic neuron i are generated by a point process 
with time dependent stochastic intensity pi(t) = g(ui(t)) that depends non- 
linearly upon the membrane potential u^. Since the membrane potential in 
turn depends on both the input and the firing history of the postsynaptic 
neuron, we write: 



We take an exponential to describe the escape across a noisy threshold, i.e, 
g(u) = po ex P f^f) where i? = —50 mV is the formal threshold, Au = 3 mV 
is the width of the threshold region and therefore tunes the stochasticity of 
the neuron, and p = 1/ms is the stochastic intensity at threshold (see Fig. HJ. 
Other choices of the escape function g are possible with no qualitative change 
of the results. For An — > 0, the model is identical to the deterministic leaky 
integrate-and-fire model with synaptic current injection |21j . 

We note that the stochastic process, defined in Eq. (JBJ) is similar to, but 
different from a Poisson process since the stochastic intensity depends on the 
set yt of the previous spikes of the postsynaptic neuron. Thus the neuron 
model has some 'memory' of previous spikes. 

2.3 Stochastic generative model 

The advantage of the probabilistic framework introduced above via the noisy 
threshold is that it is possible to describe the probability density 1 Pi(y\x) of 

1 For the sake of simplicity, we denoted the set of postsynaptic spikes from to T by y 
instead of yr- 




(4) 




(5) 
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u [mV] Iq ISI [ms] 

Figure 1: A. Escape rate g(u) = poexp f^f)- B. Firing rate of the post- 
synaptic neuron as a function of the amplitude Iq of a constant stimulation 
current (arbitrary units). C. Interspike interval (ISI) distribution for different 
input currents. 



an entire spike train 2 Y(t) = Y, t f Gv $(t — t{) (see appendix lAl for details): 

/ 

: 6XP [Jo 1 °s(/ 3 i( s l x '2/«))^( s ) ~ P»(s|x,y s )dsJ (6) 



^(2/|x) 



II pM\^y t d J ex P p».(s|x,y s )c/s 



Thus we have a generative model that allows us to describe explicitly the like- 
lihood Pi(y\x) of emitting a set of spikes y for a given input x. Moreover, since 
the likelihood in Eq. (JSJ is a smooth function of its parameters, it is straight- 
forward to differentiate it with respect to any variable. Let us differentiate 
Pi(y\x) with respect to the synaptic efficacy w^, since this is a quantity that 
we will use later on: 



glogP,(y|x) 

dwu 



o pi(s\x,y s 



[Y(s)- Pi (s\x,y s )] J2 e(s-tf)ds. (7) 



du \ u=Ui{s\x,y s )- 



where p£(s|x, y s 

In this paper, we propose three different optimal models called A, B and C 
(cf. Table [TJ. The models differ in the stimulation paradigm and the specific 
task of the neuron. The common idea behind all three approaches is the notion 
of optimal performance. Optimality is defined by an objective function L that 
is directly related to the likelihood formula of Eq. Eland that can be maximized 
by changes of the synaptic weights. Throughout the paper, this optimization 
is done by a standard technique of gradient ascent: 



Am 



1.7 



a 



dL 

Own 



(8) 



2 Capital Y is the spike train generated by the ensemble (lower case) y. 



5 



with a learning rate a. Since the three models correspond to three different 
tasks, they have a slightly different objective function. Therefore, gradient 
ascent yields slightly different strategies for synaptic update. In the following 
we start with the simplest model with the aim to illustrate the basic principles 
that generalize to the more complex models. 

3 Results 

In this section we present synaptic updates rules derived by optimizing the 
likelihood of postsynaptic spike firing at some desired firing time t des . The 
essence of the argument is introduced in a particularly simple scenario, where 
the neuron is stimulated by one presynaptic spike and the neuron is inactive 
except at the desired firing time t des . This is the raw scenario that is further 
developed in several different directions. 

Firstly, we may ask the question of how the postsynaptic spike at the de- 
sired time t des is generated: (i) The spike could simply be given by a supervisor. 
As always in maximum likelihood approaches, we then optimize the likelihood 
that this spike could have been generated by the neuron model (i.e., the genera- 
tive model) given the known input; (ii) the spike could have been generated by 
an strong current pulse of short duration applied by the supervisor (teaching 
input). In this case the a priori likelihood that the generative model fires at or 
close to the desired firing time is much higher. The two conceptual paradigms 
give slightly different results as discussed in scenario A. 

Secondly, we may, in addition to the spike at the desired time t des allow for 
other postsynaptic spikes generated spontaneously. The consequences of spon- 
taneous activity for the STDP function are discussed in scenario B. Thirdly, 
instead of imposing a single postsynaptic spike at a desired firing time t des , 
we can think of a temporal coding scheme where the postsynaptic neuron re- 
sponds to one (out of M) presynaptic spike patterns with a desired output 
spike train containing several spikes while staying inactive for the other M — 1 
presynaptic spike patterns. This corresponds to a pattern classification task 
which is the topic of scenario C. 

Moreover, optimization can be performed in an unconstrained fashion or 
under some constraint. As we will see in this section, the specific form of the 
constraint influences the results on STDP, in particular the strength of synaptic 
depression for 'pre after post' timing. To emphasize this aspect, we discuss two 
different constraints. The first constraint is motivated by the observation that 
neurons have a preferred working point defined by a typical mean firing rate 
that is stabilized by homeostatic synaptic processes |Hj . Penalizing deviations 
from a target firing rate is the constraint that we will use in scenario B. 
For very low target firing rate, the constraint reduces to the condition of 'no 
activity' which is the constraint implemented in scenario A. 

The second type of constraint is motivated by the notion of STDP itself: 
changes of synaptic plasticity should depend on the relative timing of pre- and 
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postsynaptic spike firing and not on other factors. If STDP is to be imple- 
mented by some physical or chemical mechanisms with finite time constants, we 
must require the STDP function to be local in time, i.e., the amplitude of the 
STDP function approaches zero for large time differences. This is the temporal 
locality constraint used in scenario C . While the unconstrained optimization 
problems are labeled with the subscript u (A u , B u , C u ), the constrained prob- 
lems are marked by the subscript c (A c , B c , C c ) (c.f Table HJ). 



Unconstrained scenarios 


Constrained scenarios 


A u 


Postsynaptic spike imposed 
L Au = log(p(t dcs )) 


A c 


No activity 

L Ac = L Au — Jq p{t)dt 




Postsynaptic spike imposed 
+ spontaneous activity 
L B » = log(p(t des )) 


B c 


Stabilized activity 




Postsynaptic spike 
patterns imposed 

=log(n P i (tf < |x i )JJPi(0|x fc )wb: j 

\ i k^i ) 


c c 


Temporal locality 
constraint 

L c a = L c Uj Paa , = a5AA , (- A _ f y 



Table 1: Summary of the optimality criterion L for the three unconstrained 
scenarios (A n , B u , C u ) and the three constrained scenarios (A c , B c , C c ). The 
constraint for scenario C is not included in the likelihood function L Cc itself, 
but rather in the deconvolution with a matrix P that penalizes quadratically 
the terms that are non-local in time. See appendixOfor more details. 



3.1 Scenario A: "one postsynaptic spike imposed" 

Let us start with a particularly simple model which consists of one presynaptic 
neuron and one postsynaptic neuron (c.f. Fig.EK). Let us suppose that the task 
of the postsynaptic neuron i is to fire a single spike at time t des in response to 
the input which consists of a single presynaptic spike at time t pre , i.e. the input 
is x = {t pre } and the desired output of the postsynaptic neuron is y — {t des }. 
Since there is only a single pre- and a single postsynaptic neuron involved, we 
drop in this section the indices j and i of the two neurons. 

3.1.1 Unconstrained scenario A u : 1 spike at t des 

In this subsection, we assume that the postsynaptic neuron has not been active 
in the recent past, i.e. refractory effects are negligible. In this case, we have 
jo(t|x, y t ) = p(£|x) because of the absence of previous spikes. Moreover, since 
there is only a single presynaptic spike (i.e. x = {t pre }), we write p(t|£ pre ) 
instead of p(t|x). 
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Figure 2: A. Scenario A: a single presynaptic neuron connected to a postsy- 
naptic neuron with a synapse of weight w. B. Optimal weight change given by 
Eq. (fTU|) for the scenario A u . This weight change is exactly the mirror image 
of an EPSP. 



Since the task of the postsynaptic neuron is to fire at time t des , we can 
define the optimality criterion L Au as the log-likelihood of the firing intensity 
at time t des , i.e. 

L A » = log(p(t dcs |t prc )) (9) 
The gradient ascent on this function leads to the following STDP function: 

where p'(t\t pie ) = &\ u=u nu P r:e\. Since this optimal weight change Aw Au can 
be calculated for any presynaptic firing time t pre , we get a STDP function 
which depends on the time difference At = t prc — t dcs (c.f. Fig. 133) . As we can 
see directly from Eq. (jlOj) . the shape of the potentiation is exactly a mirror 
image of an EPSP. This result is independent of the specific choice of the 
function g{u). 

The drawback of this simple model becomes apparent, if the STDP function 
given by Eq. ^TO] is iterated over several repetitions of the experiment. Ideally, 
it should converge to an optimal solution given by Aw Au = in Eo. fTUjl. 
However, the optimal solution given by Aw Au = is problematic: for At < 0, 
the optimal weight tends towards oo whereas for At > 0, there is no unique 
optimal weight (Aw Au = 0, Vw). The reason of this problem is, of course, 
that the model describes only potentiation and includes no mechanisms for 
depression. 

3.1.2 Constrained scenario A c : "no other spikes than at t des " 

In order to get some insight of where the depression could come from, let us 
consider a small modification of the previous model. In addition to the fact 
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that the neuron has to fire at time t des , let us suppose that it should not 
fire anywhere else. This condition can be implemented by an application of 
Eq. (JHJ) to the case of a single input spike x = {t pre } and a single output spike 
y = {t des }. In terms of notation we set P(y\x) = P(t des |t pre ) and similarly 
p(s\x,y) = p(s\t pre , t des ) and use Eq. (jOJ) to find: 



P(t des \t prc ) = p(t des \t pre ) exp 



( T p(s\t pie ,t des )ds 
Jo 



Note that for s < t des , the firing intensity does not depend on t des , hence 
p(s|t pre , t dcs ) = p(s\t pve ) for s < t des . We define the objective function L Ac as 
the log-likelihood of generating a single output spike at time t des , given a single 
input spike at i pre . Hence, with Eq. (jllj) : 



L = \og(P(t des \t pTC )) 

= log(p(t des \t pre )) - [ T p(s\t prc ,t dcs )ds (12) 

Jo 

and the gradient ascent Aw Ac = dL Ac /dw rule yields 

AwAc = ^^j^ ~ - jf P( s \ tPrC > ~ tPTe )ds (13) 

Since we have a single postsynaptic spike at t des , Eq. (|T3|) can directly be 
plotted as a STDP function. In Fig. |3] we distinguish two different cases. In 
Fig. we optimize the likelihood L Ac in the absence of any teaching input. 
To understand this scenario we may imagine that a postsynaptic spike has oc- 
curred spontaneously at the desired firing time t des . Applying the appropriate 
weight update calculated from Eq. (|T3*|) will make such a timing more likely 
the next time the presynaptic stimulus is repeated. The reset amplitude r} 
has only a small influence. 

In Fig. 133 we consider a case where firing of the postsynaptic spike at 
the appropriate time was made highly likely by a teaching input of duration 
AT = 1 ms centered around the desired firing t des . The form of the STDP 
function depends on the amount r] Q of the reset. If there is no reset r] = 0, 
the STDP function shows strong synaptic depression of synapses that become 
active after the postsynaptic spike. This is due to the fact that the teaching 
input causes an increase of the membrane potential that decays back to rest 
with the membrane time constant r m . Hence the window of synaptic depression 
is also exponential with the same time constant. Qualitatively the same is 
true, if we include a weak reset. The from of the depression window remains 
the same, but its amplitude is reduced. Only for strong reset to or below 
resting potential the effect inverses. A weak reset is standard in applications 
of integrate-and-fire models to in vivo data and is one of the possibilities to 
explain the high coefficient of variation of neuronal spike trains in vivo fSJ HH] • 

A further property of the STDP functions in Fig. El is a negative offset for 
|t pre — t dos | — > oo. The amplitude of the offset can be calculated for w ~ and 
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At >, i.e. Aw ~ — p'( M rest) Io° €(s)ds. This offset is due to the fact that we 
do not want spikes at other times than t des . As a result, the optimal weight 
w* (i.e. solution of Aw Au ), should be as negative as possible (w* — > — oo or 
w* — > w mm in the presence of a lower bound) for At > or At 0. 



A without teaching B with teaching 




t pre _ f des [ mg ] ^pre _ t des [ mg ] 

Figure 3: Optimal weight adaptation for scenario A c given by Eq. (JJJ) in the 
case of a teaching signal (A) and in the absence of a teaching signal (B). The 
weight change in the post-before-pre region is governed by the spike afterpo- 
tential UApif) = T)(t) + Wteach(£)- The duration of the teaching input is AT = 1 
ms. The amplitude of the current Iq is chosen so that max 4 M teac h(t) = 5 mV. 
Wrest is chosen such that the spontaneous firing rate g(u Test ) matches the desired 
firing rate 1/T, i.e. w rest = Am log S + 9 ~ —60 mV. The weight strength is 
w — 1. 



3.2 Scenario B: "spontaneous activity" 

The constraint in Scenario A c of having strictly no other postsynaptic spikes 
than the one at time t des may seem artificial. Moreover, it is this constraint 
which leads to the negative offset of the STDP function discussed at the end 
of the previous paragraph. In order to relax the constraint of "no spiking" , we 
allow in scenario B for a reasonable spontaneous activity. As above, we start 
with an unconstrained scenario B u before we turn to the constrained scenario 
B c . 

3.2.1 Unconstrained scenario B u : maximize the firing rate at t des 

Let us start with the simplest model which includes spontaneous activity. Sce- 
nario B u is the analog of the model A u , but with two differences. 

First, we include spontaneous activity in the model. Since p(i|x, y t ) de- 
pends on the spiking history for any given trial, we have to define a quantity 
which is independent of the specific realizations y of the postsynaptic spike 
train. 
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Secondly, instead of considering only one presynaptic neuron, we consider 
N = 200 presynaptic neuron, each of them emitting a single spike at time 
tj = jSt, where 5t = 1 ms (see Fig. HJ\). The input pattern will therefore be 
described by the set of delayed spikes x = {xj = {tj},j = 1, . . . , iV}. As long 
as we consider only a single spatio-temporal spike pattern in the input, it is 
always possible to relabel neurons appropriately so that neuron j + 1 fires after 
neuron j. 

A B 

0.35 1 1 1 1 1 




100 



Figure 4: Scenario B. A. N = 200 presynaptic neurons are firing one after 
the other at time tj = jSt with St = 1 ms. B. The optimal STDP function of 
scenario B u . 



Let us define the instantaneous firing rate p(t) that can be calculated by 
averaging p[t\y t ) over all realizations of postsynaptic spike trains: 

p(t\x) = (p(t\x,yt)) ytlx . (14) 

Here the notation means taking the average over all possible configura- 

tion of postsynaptic spikes up to t for a given input x. In analogy to a Poisson 
process, a specific spike train with firing times y t = {tj, tj , . . . , tf < t} is gen- 
erated with probability P(y t \x) given by Eq. (JHJ). Hence, the average O^ix °f 
Eq. (fT4*|) can be written as follows (see appendix El for numerical evaluation of 

P(t)): 

p(*I x ) = EtW •••/ p(t\x,y t )P(yt\x)dt?...tl. (15) 

b ! jo Jo 

Analogously to the model A u , we can define the quality criterion as the log- 
likelihood L Bu of firing at the desired time t des : 

L Bu = log(p(t dcs |x)). (16) 
Thus the optimal weight adaptation of synapse j is given by 

Bu _ dp(t dcs \x)/dw 
~ Pit' 



Aw Bn = UHy 1 J/UWJ (17) 
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where 9( q^ is given by 



^^ = p\t\^e(t-t J ) + Ut\^y t )^\ogP(y\^\ , (18) 
\ogP{y) is given by Eq. © and p'{t) = 



dWj vy; 10 u .y ^4- w ^ w - \ du \u=u(t), y 

Figure shows that, for our standard set of parameters, the differences 
to scenario A n are negligible. 

Figure depicts the STDP function for various values of the parameter 
An at a higher postsynaptic firing rate. We can see a small undershoot in the 
pre-bef ore-post region. The presence of this small undershoot can be under- 
stood as follows: enhancing a synapse of a presynaptic neuron that fires too 
early would induce a postsynaptic spike that arrives before the desired firing 
time and therefore, because of refractoriness, would prevent the generation of 
a spike at the desired time. The depth of this undershoot decreases with the 
stochasticity of the neuron and increases with the amplitude of the refractory 
period (if 770 = 0, there is no undershoot). In fact, correlations between pre- 
and postsynaptic firing reflect the shape of an EPSP in the high-noise regime, 
whereas they show a trough for low noise |HH ITH] . Our theory shows that the 
pre-before-post region of the optimal plasticity function is a mirror image of 
these correlations. 

3.2.2 Constrained scenario B c : "firing rate close to z/ " 

In analogy to model A c we now introduce a constraint. Instead of imposing 
strictly no spikes at times t 7^ t des , we can relax the condition and minimize 
deviations of the instantaneous firing rate p(£|x, t des ) from a reference firing 
rate vq. This can be done by introducing into Eq. (fTH|) a penalty term Pb 
given by 

/ 1 f T (p(t|x,t des ) -z/ ) 2 \ 

For small a, deviations from the reference rate yields a large penalty. For 
a — > 00, the penalty term has no influence. The optimality criterion is a 
combination of a high firing rate p at the desired time under the constraint 
of small deviations from the reference rate Uq. If we impose the penalty as a 
multiplicative factor and take as before the logarithm, we get: 

L Bc =log(p(t dcs |x)P B ) (20) 
Hence the optimal weight adaptation is given by 



Aw > = V-M " ^ h W(|x ' ( > ^ ^* |x ' 1 )dt (21) 

Since in scenario B each presynaptic neuron j fires exactly once at time tj = jSt 
and the postsynaptic neuron is trained to fire at time t Aes , we can interpret the 
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weight adaptation Awf c of Eq. (J2Tj) as a STDP function Aw Bc which depends 
on the time difference At = t pm — t des . Fig. El shows this STDP function for 
different values of the free parameter a of Eq. (|19)l . The higher the standard 
deviation a, the less effective is the penalty term. In the limit of o — > oo, the 
penalty term can be ignored and the situation is identical to that of scenario 
B u . 



A B 




Figure 5: A. The optimal STDP functions of scenario B n for different level 
of stochasticity described by the parameter Am. The standard value (Am = 3 
mV) is given by the solid line, decreased noise (An = 1 mV and Am = 0.5 
mV) are indicated by dot-dashed and dashed lines respectively . In the low- 
noise regime, enhancing a synapse which fires slightly too early can prevent 
the firing at the desired firing time t des due to refractoriness. To increase the 
firing rate at t des it is thence advantageous to decrease the firing probability 
some time before t des . Methods: For each value of Am, the initial weight wq 
are set such that the spontaneous firing rate is p = 30Hz. The amplitude of 
the STDP function is multiplied by Am in order to make them comparable. 
Reset: r] = —5 mV. B. Scenario B c . Optimal STDP function for scenario B c 
given by Eq. (f2T]) for a teaching signal of duration AT = 1 ms. The maximal 
increase of the membrane potential after 1 ms of stimulation with the teaching 
input is max t Uteach(^) = 5 mV. Synaptic efficacies Wy are initialized such that 
m = —60 mV which gives a spontaneous rate of p = u = 5 Hz. Standard 
noise level: Am = 3 mV. 



3.3 Scenario C: "pattern detection" 

3.3.1 Unconstrained scenario C u : "spike pattern imposed" 

This last scenario is a generalization of the scenario A c . Instead of restricting 
the study to a single pre- and postsynaptic neuron, we consider iV presynaptic 
neurons and M postsynaptic neurons (see Fig. EJ). The idea is to construct 
M independent detector neurons. Each detector neuron 2 = 1,..., M, should 
respond best to a specific prototype stimulus, say x', by producing a desired 
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spike train y l , but should not respond to other stimuli, i.e. y % = 0, Vx fc , k ^ i 
(see Fig. IHJ). The aim is to find a set of synaptic weights that maximize the 
probability that neuron % produces y l when x* is presented and produces no 
output when x fc , k ^ i is presented. Let the likelihood function L Cu be 

(m m \ 

IpiG/V) II ^(0|x fc )^ (22) 
i=l fc=l,fc^i J 

where Pj(y l |x l ) (c.f Eq. ©) is the probability that neuron i produces the 
spike train y % when the stimulus x* is presented. The parameter 7 characterizes 
the relative importance of the patterns that should not be learned compared 
to those that should be learned. We get 

M 

L c » = £ log^yV)) - 7 (log(P l (0|x fe ))) x (23) 

i=l 

where the notation (-) x k^ x i = jfzi Sfc^i means taking the average over all 
patterns other than x\ The optimal weight adaptation yields 

A«;g = -^-logCPiCy'lx*)) -7/^log(P(0|x fe ))\ (24) 

The learning rule of Eq. (|24|) gives the optimal weight change for each 
synapse and can be evaluated after presentation of all pre- and postsynaptic 
spike patterns, i.e. it is a "batch" update rule. Since each pre- and postsynaptic 
neuron emit many spikes in the interval [0, T], we can not directly interpret 
the result of Eq. ()24j) as a function of the time difference At = t pre — t post as 
we did in scenario A or B. 




Figure 6: Scenario C. N presynaptic neurons are fully connected to M post- 
synaptic neurons. Each postsynaptic neuron is trained to respond to a specific 
input pattern and not respond to M — 1 other patterns as described by the 
objective function of Eq. (j22J)- 



14 



pattern to be learned 
400^ 



X 200 



pattern not to be learned 
400 



X 200 



50 100 150 200 



50 100 150 200 



1000 



500 




50 100 150 200 



50 100 150 200 



^0.5 





- - — 1 











0.1 



cx0.05 



50 100 150 
time [ms] 



200 



50 100 150 
time [ms] 



200 



Figure 7: Pattern detection after learning. Top. The left raster plot represents 
the input pattern the i th neuron has to be sensitive to. Each line corresponds 
to one of the iV = 400 presynaptic neurons. Each dot represents an action 
potential. The right figure represents one of the patterns the i th neuron should 
not respond to. Middle. The left raster plot corresponds to 1000 repetitions 
of the output of neuron i when the corresponding pattern x 2 is presented. The 
right plot is the response of neuron i to one of the pattern it should not respond 
to. Bottom. The left graph represents the probability density of firing when 
pattern x l is presented. This plot can be seen as the PSTH of the middle 
graph. Arrows indicate the supervised timing neuron i learned. The right 
graph describes the probability density of firing when pattern x fc is presented. 
Note the different scales of vertical axis. 



Ideally, we would like to write the total weight change of the optimal rule 
given by Eq. (j2H) as a sum of contributions 



Awg= AW Cu (t pm -t post ), (25) 

tP re e x) 

^post £ i 



where AW Cu (t pre — t post ) is a STDP function and the summation runs over 
all pairs of pre- and postsynaptic spikes. The number of pairs of pre- and 
postsynaptic spikes with a given time shift is given by the correlation function 
which is best defined in discrete time. We assume time steps of duration 
5t = 0.5 ms. Since the correlation will depend on the presynaptic neuron j 
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t pre _ ^post [ mg ] t prc _ £ P ost [ mg ] 

Figure 8: A. Optimal weight change for scenario C u . In this case, no locality 
constraint is imposed and the result is similar to the STDP function of scenario 
A c (with r/o = and u tea ch(i) = 0) represented on Fig. El B. Optimal weight 
change for scenario C c as a function of the locality constraint characterized 
by a. The stronger the importance of the locality constraint, the narrower is 
the spike-spike interaction. For A and B: M = 20, r] Q = —5 mV. The initial 
weights Wij are chose so that the spontaneous firing rate matches the imposed 
firing rate. 



and the postsynaptic neuron i under consideration, we introduce a new index 
k = N(i — 1) + j. We define the correlation in discrete time by its matrix 
elements C^a that describe the correlation between the presynaptic spike train 
and the postsynaptic spike train Y l (t — T + A5t). For example, C 3 a = 7 
implies that 7 spike pairs of presynaptic neuron j = 3 with postsynaptic neuron 
i = 1 have a relative time shift of To — A5t. With this definition, we can rewrite 
Eq. (|23)l in vector notation (see appendix lC.il for more details): 

Aw c = CAW Cu (26) 

where Aw c = (Aw^, . . . , Awf N , Aw^, . . . , Aw^. IN ) T is the vector containing 
all the optimal weight change given by Eq. (}2"2|) and AW Cu is the vector con- 
taining the discretized STDP function with components AW^ U = AW Cu (— To+ 
A5t) for 1 < A < 2To with To = To/5. In particular, the center of th STDP 

function (i.e. t pre = tP ost ) corresponds to the index A = To. The symbol = 
expresses the fact that we want to find AW " such that Aw c is as close as 
possible to CAW Cu . By taking the pseudo- inverse C + = (C T C)~ l C T of C, 
we can invert Eq. (J2l)|) and get 

AW C " = C + Aw c (27) 

The resulting STDP function is plotted in Fig. |HK- As it was the case for the 
scenario A u , the STDP function exhibits a negative offset. In addition to the 
fact the postsynaptic neuron i should not fire at other times than the ones 
given by y l , it should also not fire whenever pattern x fc , k ^ i is presented. 
The presence of the negative offset is due to those two factors. 
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3.3.2 Constrained scenario C c : "temporal locality" 

In the previous paragraph, we obtained a STDP function with a negative offset. 
This negative offset does not seem realistic because it implies that the STDP 
function is not localized in time. In order to impose temporal locality (finite 
memory span of the learning rule) we modify Eq. (|27|) in the following way 
(see appendix IC. 21 for more details): 



where P is a diagonal matrix which penalizes non-local terms. In this paper, 
we take a quadratic suppression of terms that are non-local in time. With 
respect to a postsynaptic spike at t post , the penalty term is proportional to 
(t — t post ) 2 . In matric notation, and using our convention that the postsynaptic 
spike corresponds to A = T , we have: 



The resulting STDP functions for different values of a are plotted in Fig. [HfB. 
The higher the parameter a, the more non-local terms are penalized, the nar- 
rower is the STDP function. 

The STDP function of Fig. has been derived using a matrix inversion 
and an ad- hoc locality constraint. In order to assess th errors induced by 
our approach, we simulate an update Au^| c = J2tP ic ex z J2tp° st ey i AW Cc (t pre — 
t post ) with the STDP function AW Cc and compare this update with the exact 
optimal result Awfj from Eq. (|24|). Fig. El shows how well the reconstructed 
weight update Aw rec = CAW c describes the effective weight change Aw c 
given by Eq. (J24"j) for various values of the locality parameter a. 

Fig. ITTM - shows the STDP functions for various number of patterns M. No 
significant change can be observed for different numbers of input patterns M. 
This is due to the appropriately chosen normalization factor 1/(M — 1) in the 
exponent of Eq. (|2~2"Jl . 

The target spike trains y % have a certain number of spikes during the 
time window T, i.e. they set a target value for the mean rate. Let u post = 
Ym Jo y l (t)dt be the imposed firing rate. Let wq denote the amplitude 
of the synaptic strength such that the firing rate p WQ given by those weights 
is identical to the imposed firing rate: p m = i> post . If the actual weights are 
smaller than wq, almost all the weights should increase whereas if they are big- 
ger than Wo, depression should dominate (c.f FigllOb). Thus the exact form of 
the optimal STDP function depends on the initial weight value Wq. Alterna- 
tively, homeostatic process could assure that the mean weight value is always 
in the appropriate regime. 



AW Gc 



{C T C + P)~ 1 C T Aw° 



(28) 




(29) 
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A a = B a = 0.04 C a = 0.4 




Aw opt Aw opt Aw opt 



Figure 9: Correlation plot between the optimal synaptic weight change 
Aw opt = Aw Cu and the reconstructed weight change Aw rec = CAW c using 
the temporal locality constraint. A. No locality constraint, i.e. a = 0. Devi- 
ations from the diagonal are due to the fact that the optimal weight change 
given by Eq. (|24jl can not be perfectly accounted for the sum of pair effects. 
The mean deviations are given by Eq. (|41j) . B. A weak locality constraint 
(a = 0.04) almost does not change the quality of the weight change recon- 
struction. C. Strong locality constraint (a = 0.4). The horizontal lines arise 
since most synapses are subject to a few strong updates induced by pairs of 
pre- and postsynaptic spike times with small time shifts. 
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ipre ^post t pre ^post 

Figure 10: A. Optimal STDP function as a function of the number of input 
patterns M. (a = 0.04, N = 400) B. Optimal weight change as a function of 
the the weight w. If the weights are small (dashed line) potentiation dominates 
whereas if they are small (dotted line) depression dominates. 



4 Discussion 

Our approach is based upon the maximization of the probability of firing at 
desired times t des with or without constraints. From the point of view of 
machine learning, this is a supervised learning paradigm implemented as a 
maximum likelihood approach using the spike response model with escape noise 
as a generative model. Our work can be seen as a continuous-time extension 
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Unconstrained scenarios 


Constrained scenarios 


A u 


TiTP - hpinTP- Tin c/ 

Ul C ' UC-IUI C ■ jJUOo 

LTP ~ EPSP 


A c 


Tin Q f- flPTnTP-TiTP 
JJUo L UCJ U 1 C JJ 1 C 

LTD (or LTP) ~ spike afterpot. 


B a 


pre-bej "ore-post 

LTP/LTD ~ reverse correlation 


B c 


post-before-pre 

LTD ~ increased firing rate 


c u 


pre-bej "ore-post 
LTP ~ EPSP 

LTD ~ background patterns 


c c 


post-before-pre 
LTD ~ background patterns 
~ temporal locality 



Table 2: Main results for each scenario. 

of the maximum likelihood approach proposed by Barber j^j. In the following 
our discussion will mainly focus on the relation to neuroscience (see table |2J). 

From the point of view of neuroscience, paradigms of unsupervised or re- 
inforcement are probably much more relevant than the supervised scenario 
discussed here. For example, the gradient ascent approach, used in Scenarios 
A and B, would yield a slightly different STDP function, depending on the 
initialization of the current set of synaptic weights - and it is not obvious how 
one synapse should know about the values of all other synapses. Moreover, it 
is not clear what a target spike train imposed by a supervisor really means. 
One possible interpretation is that the target spike train is in fact not neces- 
sarily given by a supervisor, but created spontaneously by the neuron itself. 
The learning rule will then change the weights so that the same output will 
be more likely when the same stimulus is repeated. Conditioning the weight 
changes on the presence of the a reward signal, would not change the main 
results of this paper, i.e., predictions for the optimal time course of the STDP 
function. 

In scenario A, we have mainly shown that the shape of potentiation in the 
STDP function is directly related to the time course of an EPSP. Moreover, 
we have seen that the behavior of the post-before-pre region is determined by 
the spike afterpotential (see table Elfor a result summary of the three models). 
In the presence of a teaching input and firing rate constraints, a weak reset 
of the membrane potential after the spike means that the neuron effectively 
has a depolarizing spike after potential (DAP). In experiments, DAPs have 
been observed [I3IHUl2j for strong presynaptic input. Other studies, however, 
have shown that the level of depression does not depend on the postsynaptic 
membrane potential [H]. In any case, a weak reset (i.e., to a value below 
threshold rather than to the resting potential) is consistent with the findings 
of other researchers that used integrate-and-fire models to account for the high 
coefficient of variation of spike trains in vivo fT3| HH] . 

In the presence of spontaneous activity (c.f. scenario B), we have shown 
that the pre-bef ore-post regime is sensitive to the level of noise and reflects the 
correlation between pre- and postsynaptic firing. Similarly to the observations 
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made for scenario A, a constraint on the spontaneous firing rate causes the 
optimal weight change to elicit a depression of presynaptic spikes that arrive 
immediately after the postsynaptic one. In fact, the reason of the presence 
of the depression in scenario B c is directly related to the presence of a DAP 
caused by the strong teaching stimulus. In both scenarios A and B, depression 
occurs in order to compensate the increased firing probability due to the DAP. 

In scenario C, it has been shown that the best way to adapt the weights 
so that the postsynaptic neurons can detect a specific input patterns among 
others can be described as a STDP function. This task is similar to the one 
in in the sense that a neuron is designed to be sensitive to a specific input 
pattern, but different since our work does not assume any axonal delays. 

In our framework, the definition of the objective function is closely related 
to the neuronal coding. In scenario C, we postulate that neurons emit a precise 
spike train whenever the "correct" input is presented and be silent otherwise. 
This coding scheme is clearly not the most efficient one. Another possibility 
is to impose postsynaptic neurons to produce a specific but different spike 
train for each input pattern and not only for the "correct" input. Such a 
modification of the scenario does not dramatically change the results. The 
only effect is to reduce the amount of depression and increase the amount of 
potentiation. 

Theoretical approaches to neurophysiological phenomena in general, and to 
synaptic plasticity in particular, can be roughly grouped into three categories: 
first, biophysical models that aim at explaining the STDP function from prin- 
ciples of ion channel dynamics and intracelluar processes [HHl HOI HI l27j : sec- 
ond, mathematical models that start from a given STDP function and analyze 
computational principles such as intrinsic normalization of summed efficacies 
or sensitivity to correlations in the input |2Hl 03 EH HHl EOl ESI E21 1221 5 fi- 
nally, models that derive 'optimal' STDP properties for a given computational 
task EE EH El El • Optimizing the likelihood of postsynaptic firing in 
a predefined interval, as we did in this paper, is only one possibility amongst 
others of introducing concepts of optimality [3J El E] into the field of STDP. 
Chechik uses concepts from information theory [T5J, but restricts his study to 
the classification of stationary patterns. Bohte et a!., Bell et al. and Toyoizumi 
et al. El are also using concepts from information theory, but they are 
applying them to the pre- and postsynaptic spike trains. Dayan use concepts 
of optimal filter theory ^Bj , but are not interested in precise firing of the post- 
synaptic neuron. The work of Hopfield is similar to our approach in that 
it focuses on recognition of temporal input patterns, but in our approach we 
are in addition interested in triggering postsynaptic firing with precise timing. 
Hopfield and Brody emphasize the repair of disrupted synapses in a network 
that has previously acquired its function of temporal pattern detector. 

Optimality approaches, such as ours, will never be able to make strict 
predictions about the real properties of the neuron. Optimality criteria may, 
however, help to elucidate computational principles and provide insights about 



20 



potential tasks of electrophysiological phenomena such as STDP. 
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A Probability Density of a Spike Train 

The probability density of generating a spike train y t = {tj,tf, . . . , if < t} 
with the stochastic process defined by Eq. © can be expressed as follows: 

P(y t ) = P(tl...,tf)R(t\y t ) (30) 

where P{t], . . . ,tf) is the probability density of having F spikes at times 
tj, . . . ,tf and R(t\yt) = exp (— //f p(t'\S t >)dt'^J corresponds to the probability 
of having no spikes from tf to t. Since the joint probability P(tj, . . . , tf) can 
be expressed product of conditional probabilities 

P(t},...,t?) = P(tl)l[P(t{\tt\ ...,t]) (31) 

/=2 

Eq. (j3T?j) becomes 
P(y t ) = p(t}\y t} )exp^- J\(t'\y t ,)dt^j 

J] p(t{\y t f) exp (- p{t'\y v )dt^j | exp (- ^ p(t'\y t ,)dt'^ 
J] p(t{\y t{ ) exp (- jf* p(t'\y t ,)dt'^ (32) 

B Numerical Evaluation of p(t) 

Since it is impossible to numerically evaluate the instantaneous firing rate p(t) 
with the analytical expression given by Eq. ()14)) . we have to do it in a different 
way. In fact, there are two different ways to evaluate p{t). Before going into the 
details, let us first recall that from the law of large numbers, the instantaneous 
firing rate is equal to the empirical density of spikes at time t: 

(P(t\y t )) yt = (Y(t)) Y(t) (33) 
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where Y(t) = XV p $(t — t{ ) is one realization of the postsynaptic spike 
train. Thus the first and simplest method based on the r.h.s of Eq. ()33)) is to 
build a PSTH by counting spikes in small time bins [t, t + St] over, say K = 
lO'OOO repetitions of an experiment. The second, and more advanced method, 
consists in evaluating the l.h.s. of Eq. (}3*3*j) by Monte-Carlo sampling: instead of 
averaging over all possible spike trains yt, we generate K = lO'OOO spike trains 
by repetition of the same stimulus. A specific spike train y t = {t},tj, . . . ,tf < 
t} will automatically appear with appropriate probability given by Eq. (JHJ). 
The Monte-Carlo estimation p(t) of p(t) can be written as 

ffl = 4 E P(t\v?) (34) 

" m=l 

where y™ is the m th spike train generated by the stochastic process given by 
Eq. (JSJ). Since we use the analytical expression of p(t\y™), we will call Eq. (J33|) 
a semi-analytical estimation. Let us note that the semi-analytical estimation 
p(t) converges more rapidly to the true value p(t) than the empirical estimation 
based on th PSTH. 

In the limit of a Poisson process, i.e. r] Q = 0, the semi-analytical estimation 
p(t) given by Eq. (p£fj) is equal to the analytical expression of Eq. (JHJ), since 
the instantaneous firing rate p of a Poisson process is independent of the firing 
history y t = {t],tf, . . . ,tf < t} of the postsynaptic neuron. 



C Deconvolution 



C.l Deconvolution for Spike Pairs 



With a learning rule such as Eq. (|24|). we know the optimal weight change Aiw^ 
for each synapse, but we still do not know the corresponding STDP function. 

Let us first define the correlation function Cfc(r), k = N(i — 1) + j between 
the presynaptic spike train Xj(t) = Y^t^ex*. ${t ~ ^ pre ) an d the postsynaptic 
spike train Y\t) = Yt^ey* S(t - t post ): 

c k {r) = ( T X){s)Y\s + T)ds, k = l,...,NM (35) 
Jo 

where we allow a range — T < r < T , with T T. Since the sum of 
the pair based weight change AW should be equal to the total adaptation of 
weights Awk, we can write: 

[ T ° c k (s)AW(s)ds = Aw k k = 1, . . . , NM (36) 

J-T 

If we want to express Eq. (|3*3j) in a matrix form, we need to descretize time 
in small bins St and define the matrix element 

r(A+l)<St-To 

C k A = / c k (s)ds (37) 

JA5t-T 
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Now Eq. becomes 

Aw = CAW (38) 

where Aw = (Au> n , . . . , Affi lff , A«j 21 , . . . , Awmn) T is the vector contain- 
ing all the optimal weight change and AW is the vector containing the dis- 
cretized STDP function, i.e. AW A = AW(-T + A5t), for A = 1, . . . , 2T 
with T = T /5t. 

In order to solve the last matrix equation, we have to compute the inverse 
of the non-square NM x 2T matrix C, which is known as the Moore-Penrose 
inverse (or the pseudo- inverse): 

C + = {CfCy x C T (39) 

which exists only if (C T C) _1 exists. In fact, the solution given by 

AW = C+Aw (40) 

minimizes the square distance 

D = i(CAW — Aw) 2 (41) 
C.2 Temporal Locality Constraint 

If we want to impose a constraint of locality, we can add a term in the mini- 
mization process of Eq. (j4~Tj) and define the following: 

E = D + W T PAW (42) 

where P is a diagonal matrix which penalizes non-local terms. In this 
paper, we take a quadratic suppression of terms that are non-local in time: 

P AA , = a5 AA , (A-f ) 2 (43) 



To corresponds to the index of the vector AW in Eqs. (|38j) and (|42j) for 
which t pre — t post = 0. Calculating the gradient of E given by Eq. (j4*2^1 with 
respect to Aw yields 

V AW £ = C T (C AW - Aw) + PAW (44) 

By looking at the minimal value of E, i.e. VawE = 0, we have 

AW = {C T C + P)- 1 C T Aw (45) 

By setting a = 0, we recover the previous case. 
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